RStanによる反応時間の解析：Linear Ballistic Accumulator modelを用いて
専修大学人間科学部 国里愛彦
2016年12月22日
本資料は，Stan Advent Calendar 2016に向けて作成した資料です。心理学や神経科学などで用いられる反応時間について，Linear Ballistic Accumulator modelを用いて解析する方法を解説しています。本資料は， Annis et al. (2016). Bayesian inference with Stan: A tutorial on adding custom distributionsの紹介を通して，以下の２点を目的としています。
- Linear Ballistic Accumulator modelを用いた反応時間の解析を理解する
- Stanでユーザー定義関数を用いる方法を理解する
1 Linear Ballistic Accumulator modelについて
1.1 反応時間の解析とは？
心理学や神経科学などにおいて，人や動物の行動を対象とした場合、何か刺激を提示してから反応するまでにかかる時間を解析することって結構多いかと思います。こういうのを反応時間と言います。以下の図で説明すると，二重丸が出てきたらボタンを押すという実験状況があったとして，二重丸の提示からボタン押しまにかかる反応時間を解析するというものです。二重丸が出たらボタンを押すというのはあまり意味のある実験には思えないかもしれません。実際は，提示する刺激の特性を変化させたり，反応のルールを操作することで，反応時間に違いが生じるかを調べることが多いです。よく知らないですけど，ウェブ上で複数の広告を出してからどれかをクリックするまでの時間を調べたいときも，同様かなと思います。より詳しく知りたい方は，脳科学辞典の反応時間を参照ください。
1.2 従来の反応時間の解析と逐次サンプリングモデル
従来の反応時間の解析では，群間or条件間の平均値の差をt検定や分散分析で検定するというのが多かったように思います。そして，明らかになった平均値の差が，実験的な刺激の操作や群の違いを反映していると考えます。しかし，反応時間というものをよく考えると，その中にはいろいろな成分が入っています。下の図を見てください。私達が刺激が出てから反応するまでに，(1)刺激の読み込み，(2)判断，(3)ボタンを押す動作の３つがあるように思えます。多くの場合，研究者が関心のあるのは，「判断」になります。そのため，「刺激の読み込み」と「ボタンを押す動作」のような無判断時間（Non Decision Time）と判断にかかる時間（Decision Time）に分けて考えるほうが自然かと思います。このような反応時間の生成についてのモデルとして，逐次サンプリングモデルがあり，その１つにLinear Ballistic Accumulator model（無理やり訳すなら，線形弾道蓄積モデル）があります。逐次サンプリングモデルの代表的なモデルとしてdrift-diffusion modelもありますが，今回はLinear Ballistic Accumulator modelを取り上げます（Linear Ballistic Accumulator modelは，Drift-Diffusion modelをよりシンプルにしたものです）。Linear Ballistic Accumulator modelは，Brown & Heathcote(2008)によるThe simplest complete model of choice response time: Linear ballistic accumulationで提案されたモデルになります（数理モデルの導出は本論文のAppendixを参照ください）。
1.3 Linear Ballistic Accumulator modelにおけるDecision Timeの中身
さて，関心のあるDecision Timeの中身をみていきます（以降，図の中で，Non Decision Timeは左にまとめています）。Linear Ballistic Accumulator modelでは，私たちは判断を下すまでに徐々にその選択をする上でのエビデンスを蓄積し，それが閾値まで溜まったら，判断を下すと考えます（下図の上の部分）。Drift-Diffusion modelの場合は，この蓄積過程にランダムウォークが入りますが，Linear Ballistic Accumulator modelは，まさに，戦型な弾道を描いてエビデンスを蓄積していきます。
下図を用いてもう少しLinear Ballistic Accumulator modelについて説明します。Linear Ballistic Accumulator modelの前提として，反応時間は，エビデンスの蓄積が閾値(b)まで達した時になされます。その際の，エビデンスの蓄積率（傾きに近いもの）は，ドリフト率(d)と呼びます。各試行のドリフト率(d)は，平均v，標準偏差sの正規分布に従います。また，各試行の開始点(a)は，0からA(開始点の上限)の間から一様分布できまります。Decision Timeは，(b-a)/dで求めることができる。Non decision Time(τ)は，全試行で一定と考えます。ここで，aとdは，推定するパラメータではなく，v, b, A, s, τ を推定する（sは，今回は１で固定している）。今回，閾値(b)を直接計算せず，b=k+Aと考えて，相対閾値(relative threshold: k)を推定している。
1.4 Linear Ballistic Accumulator modelにおける複数選択肢の扱い
Linear Ballistic Accumulator modelでは，複数選択肢の反応時間の解析ができます（Drift-Diffusion modelの場合は，２選択肢に限定される）。まず，Linear Ballistic Accumulator modelでは，閾値までエビデンスが蓄積されたときに，反応が出力されるという前提があります。選択肢が複数ある場合，最も最初に閾値に到達した選択肢の反応が出力されます（下図の場合，選択肢１の方が，選択肢２よりも先に閾値に到達しているので選択肢１が反応として出力されます）。なお，Linear Ballistic Accumulator modelでは，Decision Timeでエビデンスを蓄積するものをaccumulatorと呼ぶ（選択肢と言っても良いが，エビデンスの蓄積過程を明示的に扱っている感じです）。
2 StanでLinear Ballistic Accumulator modelを使うためのユーザー定義関数の設定
AnnisらのBayesian inference with Stan: A tutorial on adding custom distributionsで紹介されていた関数を改変して使った。AnnisらのStanコードは，AnnisのHPにて配布されています。なお，主な改変は，(1)Annisらのコードをstan 2.10以降の記法に対応させた（<-を=に変更，_logを_lpdfに変更など），(2)Annisらのコードでは，個々の試行の対数尤度が計算されてなかったので，試行ごとの対数尤度を計算できるように修正した（後述）。
2.1 特定のaccumulatorの確率密度関数
LBAのProbability density function（確率密度関数）を計算する関数になる。引数は，t，b，A, v_PDF, sになる。入力にたいして，PDFを返す関数になる。 PDFとは，連続確率変数が取りうるある値での相対尤度を記述する関数である。この関数は Annisらのコードをそのまま使用する（書き方はStan2.11以上に合わせた）。LBAのPDFは以下の式から計算される。
\[ f_i (t) = \frac{1}{A} [-v_i Φ(\frac{b-A-tv_i}{ts})+sϕ(\frac{b-A-tv_i}{ts})+v_i Φ(\frac{b-tv_i}{ts})-sϕ(\frac{b-tv_i}{ts})] \]
real lba_pdf(real t, real b, real A, real v_pdf, real s){
//PDF of the LBA model
real b_A_tv_ts;
real b_tv_ts;
real term_1b;
real term_2b;
real term_3b;
real term_4b;
real pdf;
b_A_tv_ts = (b - A - t * v_pdf) / (t * s);
b_tv_ts = (b - t * v_pdf) / (t * s);
term_1b = v_pdf * Phi(b_A_tv_ts);
term_2b = s * exp(normal_lpdf(fabs(b_A_tv_ts)|0,1));
term_3b = v_pdf * Phi(b_tv_ts);
term_4b = s * exp(normal_lpdf(fabs(b_tv_ts)|0,1));
pdf = (1/A) * (-term_1b + term_2b + term_3b - term_4b);
return pdf;
}2.2 特定のaccumulatorの累積密度関数
LBAのCumulative Density Function（累積密度関数）を計算する関数になる。引数は，t，b，A, v_cdf, sになる。入力にたいして，CDFを返す関数になる。 PDFとは，CDFとは，ある分布がある値までに累積した確率密度の関数である。この関数は Annisらのコードをそのまま使用する（書き方はStan2.11以上に合わせた）。LBAのCDFは以下の式から計算される。
\[ F_i (t)=1+\frac{b-A-tv_i}{A} Φ(\frac{b-A-tv_i}{ts})-\frac{b-tv_i}{A}Φ\frac{b-tv_i}{ts}+\frac{ts}{A}ϕ(\frac{b-A-tv_i}{ts})-\frac{ts}{A}ϕ(\frac{b-tv_i}{ts}) \]
real lba_cdf(real t, real b, real A, real v_cdf, real s){
//CDF of the LBA model
real b_A_tv;
real b_tv;
real ts;
real term_1a;
real term_2a;
real term_3a;
real term_4a;
real cdf;
b_A_tv = b - A - t * v_cdf;
b_tv = b - t * v_cdf;
ts = t * s;
term_1a = b_A_tv / A * Phi(b_A_tv / ts);
term_2a = b_tv / A * Phi(b_tv / ts);
term_3a = ts / A * exp(normal_lpdf(fabs(b_A_tv / ts)|0,1));
term_4a = ts / A * exp(normal_lpdf(fabs(b_tv / ts)|0,1));
cdf = 1 + term_1a - term_2a + term_3a - term_4a;
return cdf;
}2.3 特定のaccumulatorが最初に閾値に到達する累積密度関数
上で定義した，lba_cdfとlba_pdfは，時間tにおいて，i番目のaccumulatorが閾値に到達した確率密度や累積密度を出している。それを用いて，複数のaccumulatorがある状況において，どのaccumulatorが最初に閾値につくのかに興味がある（つまり，i番目のaccumulatorにおける反応時間のdefective分布に関心がある）。defective確率密度関数は，lab_cdとlba_pdfを使って以下のような式で計算される。以下の場合，i番目のaccumulatorが一番最初に閾値に到達する確率密度関数は，accumulator iの確率密度関数と，(1-それ以外のaccumulatorの累積密度関数)の総積をかけたものになる。つまり，選択肢iを時間tに選ぶ確率＊他の選択肢の累積確率の残りの総積をしている。これは，選択肢iで反応が起こる確率だけでなく，その他の選択肢で反応が起こらない確率も考慮している。
\[ PDF_i (t)=f_i (t)\prod_{j\neq i}(1-F_j (t)) \]
2.4 LBAの対数尤度関数(オリジナル)
Annisらによるlba_logは，パラメータと反応時間を用いて，どちらを選択するかをPDFで計算したうえで，その選択確率の尤度を計算する関数になる。なお，閾値(b)は，開始点の上限＋ｋ（相対閾値）であり，t(decision time)は，反応時間-non decision timeである。引数はmatrix RT, real k, real A, vector v, real s, real psiになる。この関数は，RTを行列にしており，１つの実験のブロック単位での対数尤度を計算する。これは，WAICを計算する上で，少々問題になる（WAICの計算では，各従属変数ごとの対数尤度が必要になる）。以下で，Annisらのコードを改変する。
real lba_log(matrix RT, real k, real A, vector v, real s, real psi){
real t;
real b;
real cdf;
real pdf;
vector[cols(RT)] prob;
real out;
real prob_neg;
b <- A + k;
for (i in 1: cols(RT)){
t <- RT[1, i] - psi;
if(t > 0){
cdf <- 1;
for(j in 1: num_elements(v)){
if(RT[2, i] == j){
pdf <- lba_pdf(t, b, A, v[j], s);
}else{
cdf <- lba_cdf(t, b, A, v[j], s) * cdf;
}
}
prob_neg <- 1;
for(j in 1: num_elements(v)){
prob_neg <- Phi(-v[j]/s) * prob_neg;
}
prob[i] <- pdf * (1-cdf);
prob[i] <- prob[i] / (1-prob_neg);
if(prob[i] < 1e-10){
prob[i] <- 1e-10;
}
}else{
prob[i] <- 1e-10;
}
}
out <- sum(log(prob));
return out;
}2.5 LBAの対数尤度関数（改変版）
上記の問題点から，各試行ごとに対数尤度を計算する関数を作成した。また，rstan2.11以上に対応するために，<-を=に変更したり，_logを_lpdfに変更した。なお，rstan2.11ではユーザ定義関数を生成量ブロックで使う際にバグが生じて動かなかったが，rstan2.12では改善した。最初に，閾値(b)は，開始点の上限＋ｋ（相対閾値）であり，t(decision time)は，反応時間-non decision timeとしている。Annisらのコードでは，データの引数は，RTだけだったが，このコードでは，RT（１ブロック内の反応時間と反応の３次元配列）をrt（特定の試行の反応時間）とres（特定の試行の反応）に分けている。選択した方の選択肢のPDFと選択してない方のCDF（厳密には，選択してない累積確率）を計算する。それを，prob_negでは，vがマイナスになった場合の対処をとっているらしい。prob = pdf * cdfがPDFの計算になる。最後に，out = log(prob)で対数尤度を計算している。
※ Phiは，cumulative unit normal distribution function。
real lba_lpdf(real rt, real res, real k, real A, vector v, real s, real tau){
real t;
real b;
real cdf;
real pdf;
real prob;
real out;
real prob_neg;
b = A + k;
t = rt - tau;
if(t > 0){
cdf = 1;
for(j in 1 : num_elements(v)){
if(res == j){
pdf = lba_pdf(t, b, A, v[j], s);
}else{
cdf = (1 - lba_cdf(t, b, A, v[j], s)) * cdf;
}
}
prob_neg = 1;
for(j in 1:num_elements(v)){
prob_neg = Phi(-v[j] / s) * prob_neg;
}
prob = pdf * cdf;
prob = prob / (1 - prob_neg);
if(prob < 1e-10){
prob = 1e-10;
}
}else{
prob = 1e-10;
}
out = log(prob);
return out;
}2.6 LBAのrng関数
rngは，乱数発生関数になる。パラメータを引数として，そこから反応時間と反応の乱数を作成する。出力はpredであり，pred[1]が反応時間，pred[2]が反応になる。この関数はAnnisらのコードをそのまま使用する（書き方はStan2.11以上に合わせた）。
vector lba_rng(real k, real A, vector v, real s, real psi){
int get_pos_drift;
int no_pos_drift;
int get_first_pos;
vector[num_elements(v)] drift;
int max_iter;
int iter;
real start[num_elements(v)];
real ttf[num_elements(v)];
int resp[num_elements(v)];
real rt;
vector[2] pred;
real b;
//try to get a positive drift rate
get_pos_drift = 1;
no_pos_drift = 0;
max_iter = 1000;
iter = 0;
while(get_pos_drift){
for(j in 1 : num_elements(v)){
drift[j] = normal_rng(v[j], s);
if(drift[j] > 0){
get_pos_drift = 0;
}
}
iter = iter + 1;
if(iter > max_iter){
get_pos_drift = 0;
no_pos_drift = 1;
}
}
//if both drift rates are <= 0
//return an infinite response time
if(no_pos_drift){
pred[1] = -1;
pred[2] = -1;
}else{
b = A + k;
for(i in 1 : num_elements(v)){
//start time of each accumulator
start[i] = uniform_rng(0, A);
//finish times
ttf[i] = (b-start[i]) / drift[i];
}
//rt is the fastest accumulator finish time
//if one is negative get the positive drift
resp = sort_indices_asc(ttf);
ttf = sort_asc(ttf);
get_first_pos = 1;
iter = 1;
while(get_first_pos){
if(ttf[iter] > 0){
pred[1] = ttf[iter];
pred[2] = resp[iter];
get_first_pos = 0;
}
iter = iter + 1;
}
}
return pred;
}3 RStanによるLinear Ballistic Accumulator modelを用いた解析
3.1 使用パッケージ
ワークスペースのクリアをして，以下のパッケージを読み込みます。
rm(list=ls(all=TRUE))
#グラフィカルモデル用
library(DiagrammeR)
#グラフ作成用
library(ggplot2)
library(plotly)
library(bayesplot)
#Stan用
library(rstan)
library(loo)
#LBA用データの生成
library(rtdists)3.2 グラフィカルモデル
\[ k \sim Normal(.5,1)T[0,]\\ A \sim Normal(.5,1)T[0,]\\ \tau \sim Normal(.5,.5)T[0,]\\ v[1],v[2] \sim Normal(2,1)T[0,]\\ S = 1 \hspace{30pt}※Sは定数\\ RT[t] \sim LBA(k,A,\tau,v1,s) \hspace{30pt}※選択肢１の時\\ RT[t] \sim LBA(k,A,\tau,v2,s) \hspace{30pt}※選択肢2の時 \]3.3 データの作成
#データセットの作成（t0 = non decision time）
set.seed(123)
data <- rLBA(500, A=0.5, b=1, t0 = 0.5, mean_v=c(2, 1.5), sd_v=c(1,1))## Results based on 2 accumulators/drift rates.trialLength = length(data$rt)
stanData <- list(rt=data$rt,res=data$response,LENGTH=trialLength,NUM_CHOICES=2)3.4 データの確認
table(data$response)##
## 1 2
## 316 184# 選択肢1の反応時間
data1 <- subset(data,data$response==1)
p1 <- ggplot(data1,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 70))
# 選択肢2の反応時間
data2 <- subset(data,data$response==2)
p2 <- ggplot(data2,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 70))
p <- subplot(p1, p2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.ggplotly(p)3.5 Stanコード
※なお，以降は，kazutan氏のStan Advent 2016の記事を参考に，R MarkdownにStanコード直書きしてコンパイルする方法をとっています（私のRの利用の仕方からすると，これ便利です）。
functions{
real lba_pdf(real t, real b, real A, real v, real s){
//PDF of the LBA model
real b_A_tv_ts;
real b_tv_ts;
real term_1;
real term_2;
real term_3;
real term_4;
real pdf;
b_A_tv_ts = (b - A - t*v)/(t*s);
b_tv_ts = (b - t*v)/(t*s);
term_1 = v*Phi(b_A_tv_ts);
term_2 = s*exp(normal_lpdf(b_A_tv_ts|0,1));
term_3 = v*Phi(b_tv_ts);
term_4 = s*exp(normal_lpdf(b_tv_ts|0,1));
pdf = (1/A)*(-term_1 + term_2 + term_3 - term_4);
return pdf;
}
real lba_cdf(real t, real b, real A, real v, real s){
//CDF of the LBA model
real b_A_tv;
real b_tv;
real ts;
real term_1;
real term_2;
real term_3;
real term_4;
real cdf;
b_A_tv = b - A - t*v;
b_tv = b - t*v;
ts = t*s;
term_1 = b_A_tv/A * Phi(b_A_tv/ts);
term_2 = b_tv/A * Phi(b_tv/ts);
term_3 = ts/A * exp(normal_lpdf(b_A_tv/ts|0,1));
term_4 = ts/A * exp(normal_lpdf(b_tv/ts|0,1));
cdf = 1 + term_1 - term_2 + term_3 - term_4;
return cdf;
}
real lba_lpdf(real rt, real res, real k, real A, vector v, real s, real tau){
real t;
real b;
real cdf;
real pdf;
real prob;
real out;
real prob_neg;
b = A + k;
t = rt - tau;
if(t > 0){
cdf = 1;
for(j in 1 : num_elements(v)){
if(res == j){
pdf = lba_pdf(t, b, A, v[j], s);
}else{
cdf = (1 - lba_cdf(t, b, A, v[j], s)) * cdf;
}
}
prob_neg = 1;
for(j in 1:num_elements(v)){
prob_neg = Phi(-v[j] / s) * prob_neg;
}
prob = pdf * cdf;
prob = prob / (1 - prob_neg);
if(prob < 1e-10){
prob = 1e-10;
}
}else{
prob = 1e-10;
}
out = log(prob);
return out;
}
vector lba_rng(real k, real A, vector v, real s, real tau){
int get_pos_drift;
int no_pos_drift;
int get_first_pos;
vector[num_elements(v)] drift;
int max_iter;
int iter;
real start[num_elements(v)];
real ttf[num_elements(v)];
int resp[num_elements(v)];
real rt;
vector[2] pred;
real b;
//try to get a positive drift rate
get_pos_drift = 1;
no_pos_drift = 0;
max_iter = 1000;
iter = 0;
while(get_pos_drift){
for(j in 1:num_elements(v)){
drift[j] = normal_rng(v[j],s);
if(drift[j] > 0){
get_pos_drift = 0;
}
}
iter = iter + 1;
if(iter > max_iter){
get_pos_drift = 0;
no_pos_drift = 1;
}
}
//if both drift rates are <= 0
//return an infinite response time
if(no_pos_drift){
pred[1] = -1;
pred[2] = -1;
}else{
b = A + k;
for(i in 1:num_elements(v)){
//start time of each accumulator
start[i] = uniform_rng(0,A);
//finish times
ttf[i] = (b-start[i])/drift[i];
}
//rt is the fastest accumulator finish time
//if one is negative get the positive drift
resp = sort_indices_asc(ttf);
ttf = sort_asc(ttf);
get_first_pos = 1;
iter = 1;
while(get_first_pos){
if(ttf[iter] > 0){
pred[1] = ttf[iter] + tau;
pred[2] = resp[iter];
get_first_pos = 0;
}
iter = iter + 1;
}
}
return pred;
}
}
data{
int LENGTH;
int NUM_CHOICES;
vector[LENGTH] rt;
vector[LENGTH] res;
}
parameters {
real<lower=0> k;
real<lower=0> A;
real<lower=0> tau;
vector<lower=0>[NUM_CHOICES] v;
}
transformed parameters {
real s;
s = 1;
}
model {
k ~ normal(.5,1)T[0,];
A ~ normal(.5,1)T[0,];
tau ~ normal(.5,.5)T[0,];
for(n in 1:NUM_CHOICES){
v[n] ~ normal(2,1)T[0,];
}
for(m in 1:LENGTH){
rt[m] ~ lba(res[m],k,A,v,s,tau);
}
}
generated quantities {
vector[2] pred;
vector[LENGTH] log_lik;
pred = lba_rng(k,A,v,s,tau);
for(i in 1 : LENGTH){
log_lik[i] = lba_lpdf(rt[i] | res[i], k, A, v, s, tau);
}
}3.6 Rstanによるサンプリング
# オプション設定
rstan_options(auto_write=TRUE)
options(mc.cores = parallel::detectCores())
hmcIter = 1000
hmcChains = 4
hmcWarmup = 300
hmcThin = 2
# サンプリング
fit <- sampling(lbaModel,
seed = 1234,
data = stanData,
warmup = hmcWarmup,
iter = hmcIter,
chains = hmcChains,
thin =hmcThin)3.7 結果
3.7.1 収束判定
まあ，収束してそうですね。収束判定には，bayesplotを活用しています。このページが参考になりました。 なお，定数のSは除外しています。
3.7.1.1 トレースプロット
posterior1 <- extract(fit, inc_warmup = TRUE, permuted = FALSE)
color_scheme_set("mix-blue-pink")
p <- mcmc_trace(posterior1, pars = c("k", "A", "tau", "v[1]", "v[2]"), n_warmup = 125,
facet_args = list(nrow = 2, labeller = label_parsed))
p + facet_text(size = 15)
plot(p)
3.7.1.2 \(\hat{R}\)の確認
color_scheme_set("brightblue")
rhats <- rhat(fit,par = c("k", "A", "tau", "v[1]", "v[2]"))
print(rhats)## k A tau v[1] v[2]
## 1.002922 1.000405 1.002499 1.001867 1.002770mcmc_rhat(rhats) + yaxis_text()
3.7.1.3 有効サンプルサイズ
neffRatios <- neff_ratio(fit,par = c("k", "A", "tau", "v[1]", "v[2]"))
print(neffRatios)## k A tau v[1] v[2]
## 0.3043613 0.3113776 0.3297243 0.4265390 0.3862560mcmc_neff(neffRatios) + yaxis_text()
3.7.1.4 自己相関
posterior2 <- as.matrix(fit)
mcmc_acf(posterior2, pars = "k", lags = 10)
mcmc_acf(posterior2, pars = "A", lags = 10)
mcmc_acf(posterior2, pars = "tau", lags = 10)
mcmc_acf(posterior2, pars = "v[1]", lags = 10)
mcmc_acf(posterior2, pars = "v[2]", lags = 10)
3.7.2 No-U-Turn Samplerの診断
3.7.2.1 Divergent transitions
mcmc_nuts_divergence(nuts_params(fit), log_posterior(fit))
3.7.2.2 Energy and Bayesian fraction of missing information
mcmc_nuts_energy(nuts_params(fit),merge_chains = FALSE)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
3.7.3 推定結果
extractFit <- rstan::extract(fit)
# MAP推定値計算関数
mapEstimate <- function(z){
density(z)$x[which.max(density(z)$y)]
}
kMap <- mapEstimate(extractFit$k)
AMap <- mapEstimate(extractFit$A)
tauMap <- mapEstimate(extractFit$tau)
v1Map <- mapEstimate(extractFit$v[,1])
v2Map <- mapEstimate(extractFit$v[,2])
print(fit,par = c("k", "A", "tau", "v[1]", "v[2]","s"))## Inference for Stan model: a3ca8da18d6e18727882a51b5a14187a.
## 4 chains, each with iter=1000; warmup=300; thin=2;
## post-warmup draws per chain=350, total post-warmup draws=1400.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## k 0.85 0.01 0.24 0.45 0.66 0.84 1.01 1.34 426 1
## A 0.39 0.01 0.22 0.02 0.20 0.39 0.56 0.78 436 1
## tau 0.44 0.00 0.04 0.36 0.41 0.44 0.47 0.51 462 1
## v[1] 2.46 0.01 0.18 2.11 2.33 2.44 2.57 2.84 597 1
## v[2] 1.98 0.01 0.19 1.62 1.85 1.97 2.10 2.41 541 1
## s 1.00 0.00 0.00 1.00 1.00 1.00 1.00 1.00 1400 NaN
##
## Samples were drawn using NUTS(diag_e) at Wed Dec 14 18:55:27 2016.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).print("MAP推定値")## [1] "MAP推定値"print(paste("k=",kMap,",A=",AMap,",tau=",tauMap,",v1=",v1Map,",v2=",v2Map))## [1] "k= 0.845690572305835 ,A= 0.420898353864592 ,tau= 0.429237974027454 ,v1= 2.41385166405983 ,v2= 1.9288280320449"mcmc_areas(posterior2, pars = c("k", "A", "tau", "v[1]", "v[2]"), prob = 0.95)
3.7.4 事後予測チェック
可能なら，サブプロットごとにタイトルをつけたいです。
pred_rt = posterior2[,'pred[1]']
pred_resp = posterior2[,'pred[2]']
# 500サンプルを抽出
postPredData <- data.frame(pred_rt[1:500],pred_resp[1:500])
names(postPredData) <- c("predRt","predRes")
table(data$response)##
## 1 2
## 316 184table(postPredData$predRes)##
## 1 2
## 304 196# データの選択肢1の反応時間
data1 <- subset(data,data$response==1)
p1 <- ggplot(data1,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 70))
# データの選択肢2の反応時間
data2 <- subset(data,data$response==2)
p2 <- ggplot(data2,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 70))
# 事後予測データの選択肢1の反応時間
postPredData1 <- subset(postPredData,postPredData$predRes==1)
p3 <- ggplot(postPredData1,aes(x = predRt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 70))
# 事後予測データの選択肢2の反応時間
postPredData2 <- subset(postPredData,postPredData$predRes==2)
p4 <- ggplot(postPredData2,aes(x = predRt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 70))
p <- subplot(p1,p2,p3,p4)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.ggplotly(p)3.7.5 情報量基準(WAICとLOOCV)の算出
@berobero11氏のStan Advent 2016を参考にしつつWAICとLOOCVの算出をしてみた。今回は，あまり意味ないけど，今後モデル比較する場合には，活用できる。
loo::loo(extractFit$log_lik)$looic/(2*trialLength)## [1] 0.02337576loo::waic(extractFit$log_lik)$waic/(2*trialLength)## [1] 0.02336874